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We consider the time-reversal-invariant Hofstadter-Hubbard model which can be realized in cold 
atom experiments. In these experiments, an additional staggered potential and an artificial Rashba- 
type spin-orbit coupling are available. Without interactions, the system exhibits various phases such 
as topological and normal insulator, metal as well as semi-metal phases with two or even more Dirac 
cones. Using a combination of real- space dynamical mean- field theory and analytical techniques, 
we discuss the effect of on-site interactions and determine the corresponding phase diagram. In 
particular, we investigate the semi-metal to antiferromagnetic insulator transition and the stability 
of different topological insulator phases in the presence of strong interactions. We compute spectral 
functions which allow us to study the edge states of the strongly correlated topological phases. 



PACS numbers: 67.85-d, 37.10.Jk 

Introduction. — Ultracold quantum gases trapped in 
optical lattice potentials provide insight into strongly 
correlated condensed matter systems. Examples are the 
Mott-insulator-superfluid transition, the dynamics of the 
Hubbard model after a quench of parameters, and the 
simulation of quantum magnetism [1] . Striking is the pre- 
cise experimental control over almost all system parame- 
ters, including the particle-particle interaction strength. 
Simulating more traditional electronic condensed matter 
systems, however, is complicated by the fact that cold 
atoms are charge neutral, and their center of mass mo- 
tion is thus not affected by external magnetic or electric 
fields (apart from trapping potentials). An experimental 
breakthrough was thus the engineering of so-called "arti- 
ficial" gauge fields, which give rise to effective magnetic 
or electric fields for the neutral particles [2]. Remark- 
ably, they may even be generalized to simulate spin-orbit 
couplings or non-Abelian fields [3]. The effective electro- 
magnetic fields and couplings can be large, which allows, 
for example, realization of the quantum (spin) Hall effect 
in a completely new experimental context [4-8]. 

The underlying idea of realizing time-reversal-invariant 
two-dimensional (2D) topological phases with cold atoms 
is as simple as it is fundamental [5, 7]. Consider the 
(integer) quantum Hall effect (QHE) on a 2D square 
lattice where an external magnetic field along the z di- 
rection breaks time-reversal and translational symmetry. 
The single particle spectrum for arbitrary magnetic field 
strength - having the shape of a butterfiy - was first 
computed by Douglas Hofstadter [9] since then referred 
to as the Hofstadter butterfiy. If the magnetic fiux per 
plaquette is a rational number a = p/q^ in units of the 
Dirac fiux quantum $o = h/e^ the system remains trans- 
lationally invariant with an enlarged unit cell of q lattice 
sites. The spectrum consists of q energy bands and in 
all energy gaps one finds a finite Chern number C and 



correspondingly \C\ chiral edge modes per edge. Inter- 
estingly, for even values of q the system is a semi-metal 
at half-filling and exhibits q Dirac cones. 

To restore time-reversal symmetry we can imagine 
applying a magnetic field in the z direction that only 
couples to the up-spins and a second field of the same 
strength but opposite direction that only couples to the 
down-spins. We thus end up with a spinful and time- 
reversal- invariant (TRI) version of the fundamental Hof- 
stadter problem. Remarkably, such a scenario is feasible 
using cold atoms in artificial gauge fields [7, 8]. Thus, the 
semi-metallic Dirac dispersion for even q becomes a gen- 
eralization of graphene with a tunable number of Dirac 
cones. Energy gaps which were crossed by a single chi- 
ral edge mode in the QHE setup are now traversed by a 
helical Kramer's pair of edge states, corresponding to a 
topological insulator phase. Note that one can use the 
same Gedanken experiment to construct the Kane-Mele 
model [4] from two time-reversed copies of Haldane's hon- 
eycomb model [10]. The Kane-Mele model with addi- 
tional Hubbard interaction has recently been intensively 
studied [11], in contrast to the Hofstadter problem. 

In this Letter we study the effect of interactions in 
the (TR- invariant) Hofstadter-Hubbard model using real- 
space dynamical mean-field theory (RDMFT) [12]. We 
explain our numerical results using analytical arguments. 
We consider interaction effects on both (semi-) metallic 
and gapped topological phases. Although Z2 topolog- 
ical insulators are known to be robust against disor- 
der [13, 14], rigorous and general results about the fate 
of topological insulators in the presence of Coulomb 
or Hubbard interactions are limited [15]. Some three- 
dimensional materials of the iridate family are possi- 
ble candidates for systems where strong spin-orbit cou- 
pling and Coulomb interactions compete [16, 17]. In 
two dimensions, however, topological insulator phases 
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have so far only been found in HgTe/CdTe quantum 
wells [18, 19], where Coulomb interactions seem to be 
negligible. 

Interacting TRI Hofstadter problem. — The TRI 
Hofstadter-Hubbard model is described by the Hamil- 
tonian 



(1) 



where = (c^^, c^^) at lattice site j = (x, ?/), is a Pauli 
matrix and x = (1,0), y = (0, 1) are unit vectors, tx {ty) 
is the hopping amplitude in x{y-) direction. We focus 
on isotropic hopping tx = ty = t here, and express all 
energies in units of t = 1. The value of a determines the 
strength of the (artificial) magnetic field for either spin 
species which penetrates a lattice plaquette in units of 
the Dirac fiux quantum. The on-site interaction strength 
U can be experimentally tuned by Feshbach resonances 
and by adjusting the lattice depth. For U = this model 
was studied in Ref. 7 (for experimental details see Sup- 
plemental Material). 

We first consider the TRI Hofstadter-Hubbard prob- 
lem for general a = p/g at half filling. For q odd the 
system is metallic with a nested Fermi surface, and anti- 
ferromagnetic Neel order occurs for infinitesimally small 
interaction U = 0+ as for the ordinary square lattice. For 
q even the situation is very different because the system is 
a semimetal (SM) at half filling [20]. The noninteracting 
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FIG. 1. (color online). Magnetization m = — n^. in the 
Neel state is plotted versus interaction strength U. We show 
results for a = 1/2 (blue), 1/4 (red), 1/6 (magenta), 1/8 
(black), and 1/10 (green) (curves from left to right). Inset: 
Fermi velocity 27tvf (red symbols) for different a = 1/q is 
shown versus q. Uc (blue symbols) obtained within RDMFT 
versus q is also shown. Magenta (upper) line is a fit of t'i? to 
oc 1/q and cyan (lower) line ol Uc to oc 1/q^ . Note that odd q 
denominators exhibit Uc = 0^ due to a nested Fermi surface. 



band structure exhibits q Dirac cones (with a multiplicity 
of 2 due to spin) which are separated by momentum 27r/g 
in momentum space. The a = 1/2 case is thus very sim- 
ilar to graphene (but note that the coordination number 
is 2: = 4 rather than 2: = 3). For smaller a on the other 
hand the system embodies a generalization of graphene 
with a tunable number of valleys. 

We investigate the SM-insulator transition for various 
a = 1/q {q even) within RDMFT. In Fig. 1, the magneti- 
zation is shown as a function of interaction U. The insu- 
lating phase for U > Uc antiferromagnetically (AFM) 
ordered with a magnetization pointing in the 2;-direction 
and an ordering wave vector Q = (7r,7r). We find that 
the critical value of Uc to enter the insulating and mag- 
netically ordered phase decreases for increasing q. This 
is expected from the increased scattering that can take 
place between the cones. At Uc we also observe a simul- 
taneous opening of the single particle gap. Within our 
approach we thus find no sign of an intermediate non- 
magnetic gapped phase. 

To understand the behavior of Uc{q) we make use of 
Herbut's argument [21]. Herbut considers graphene and 
studies the SM-insulator transition within a large- A/" ap- 
proach, and finds that Uc depends on 2A/', the number 
of Dirac cones {N refers to the spin degeneracy), and 
the Fermi velocity vp as Uc ~ vf/2N. As shown in de- 
tail in the Supplemental Material, we are able to match 
our results with Herbut's analysis by replacing the Fermi 
velocities and 2N = qN . In fact, from the bandstruc- 
ture at /7 = we find vp oc 1/q. Consequently, setting 
A^ = 2 for spin- 1/2 particles, Uc should exhibit a 1/g^ 
behavior which agrees very well with the RDMFT data; 
see inset of Fig. 1. We further note that we find that 
Uc{oc = 3/8) < Uc{(^ = 1/8), which is in agreement with 
the analysis above since vf{ol = 3/8) < vf{oc = 1/8). 

Specific cold atom setup. — We now consider two addi- 
tional terms in the Hamiltonian that are available in the 
cold-atom setup [7, 8]: a staggering of the optical lattice 
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FIG. 2. (color online). Real space magnetization profile m(x) 
in S^-S^ plane for a = 1/6, [/ = 5, = 0, and 7 = 0.125 
(a) and 7 = 0.25 (b), respectively. 
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FIG. 3. (color online), (a) ^F-A^c-phase diagram at 7 = [/ = 0. (b) [/-A^-phase diagram at np = 2/3 and 7 = 0. (c) 
Acc-7-pliase diagram at half filling = 1 and [7 = 0. (d) U-Xx-phase diagrams at Ep = for various values of 7 [indicated 
by arrows in part (c)] and inverse temperature /S = 20. We find (semi)metal (red), normal insulator (NI) (green, in upper left), 
topological insulator (QSH) (blue), and magnetically ordered phases (purple, in lower right). 



potential along the x direction 

Hx = ^(-1)-A,4c,- , 



(2) 



and a Rashba-like spin-orbit coupling that breaks axial 
spin symmetry. It is introduced via replacing in Eq. (15) 



tx tx exp(— z27r7cr^) . 



(3) 



We first study the effect of finite A^^ and 7 on the magnetic 
ordering. 

Tunable magnetism. — For 7 = 0, increasing A^^ in- 
creases Uc but does not change the type of magnetic or- 
der. Finite 7 does, however, change the type of magnetic 
order in general. To demonstrate this we consider fixed 
[/ = 5ata = l/6 and calculate the magnetization in 
real space for various values of 7 G [0,0.25]. The spin 
operators are defined in terms of the fermionic operators 
as usual as Sj = \c^-crcj with cr = (a^, a^, a^). We show 
the magnetization pattern for 7 = 0.125 and 7 = 0.25 
in Fig. 2 obtained within RDMFT. We obtain similar re- 
sults for other values of a and 7. For 7 = 0.125, the 
magnetization lies in the S^-S^ plane, has a periodicity 
of six (two) lattice sites along x{y) and reads explicitly 
m{x) = 5'tot {x) cos TTT/ (0, — cos( ^ + T]) , sin( ^ + 77)) with 
T] = 0.39 sin ^ COS The magnetization makes an- 
gles of {0°, 70°, 110°, 180°, 250°, 290°} in the Sy-S^ plane 
(spiral order). For 7 = 0.25, the magnetic order is given 
by m{x) = 5'tot(^)(0, 0, cosTT?/) (collinear order). Quan- 
tum fluctuations reduce the size of the magnetization 
-^tot < 1/2, which depends not only on the parameters 
a, 7 and U/t^ but is also spatially staggered for interme- 
diate values oiU/t (see Fig. 2). The staggering decreases 
for larger values of U/t. More importantly, tuning the 
parameter 7 we pass from Neel to spiral to collinear or- 
der crossing two magnetic quantum phase transitions. 

We can qualitatively understand this type of magnetic 
order by rigorously deriving a quantum spin Hamiltonian 



for even stronger interactions when charge fluctuations 
freeze out at half filling (see Supplemental Material for 
details) 
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where Ji = Atl/U . The first part describes spin exchange 
in X direction. For 7 = n/2 with n G Z we obtain a simple 
antiferromagnetic Heisenberg term. Other values of 7, 
however, break SU(2) symmetry and cause anisotropy of 
the XXZ type with 5*^ as anisotropy direction in spin 
space. For 7 7^ n/4 there is an additional Dzyaloshinskii- 
Moriya (DM) interaction term in the Y Z plane, which 
is responsible for the spiral spin order in Fig. 2(a). Spin 
exchange in the y direction is periodic with an extended 
unit cell in the x direction depending on the flux a = p/q'. 
for odd q the unit cell contains q lattice sites, but for even 
q it only contains q/2 lattice sites, reflecting second order 
perturbation theory. For instance, one finds for the tt- 
flux lattice {a = 1/2) an ordinary Heisenberg exchange 
term. For other values of a the XY term exhibits a 
modulation of its amplitude depending on a, while the Z 
term always favors AFM Ising order. The rich magnetic 
order predicted by the spin Hamiltonian is in agreement 
with our RDMFT findings. 

Topological insulators. — Let us now turn to the study 
of interaction effects on the gapped phases. For = 0, 
we distinguish the normal (NI) and topological (TI) in- 
sulating phases by calculating the Z2 invariant v using 
Hatsugai's method [22]. For U > 0, we identify the 
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phases by computing the spectral function in a cyhn- 
drical geometry using RDMFT and counting the number 
of gapless hehcal edge states crossing the bulk gap (for 
technical details we refer to the Supplemental Material). 
The TI phase exhibits an odd number of helical Kramer's 
pairs per edge while the NI phase an even number (in- 
cluding zero). Edge states are also crucial for detection 
of topological phases in cold-atom experiments, and we 
numerically study how robust they are with respect to 
interactions. In the following, we focus on fixed a = 1/6, 
which qualitatively captures all phenomena that occur in 
this system for general a = p/q. 

In the axial symmetric case of 7 = there exist TI 
phases only away from half filling, since the system is a 
(semi)metal for np = 1 (and not too large Xx^U). This 
is shown in Fig. 3(a), and is expected as the spinless Hof- 
stadter problem at a = 1/6 exhibits a QHE with Chern 
number C = ±2 for in the two energy gaps closest 
to zero and a QHE with C = ±1 for Ep in the other 
gaps. The Chern number corresponds to the number of 
chiral edge modes in an open geometry. In the present 
time-reversal invariant system we thus find an according 
number of helical Kramer's pairs within the gaps. 

For a filling of rip = 1/3, 5/3 the system is thus a TI. 
We observe this topological phase to be stable even for 
large interactions up to /7 = 10. We can induce a NI- 
TI phase transition in the other gap for np = 2/3,4/3 
by applying a large enough staggered lattice potential 
> 1 [see Fig. 3(a)]. Fixing np = 2/3 we now turn on 
interactions, and observe that this phase is quite stable 
as shown in Fig. 3(b). Eventually, large enough inter- 
actions reverse the effect of the staggering potential and 
drive the system into the NI phase. Note that a static 
Hartree-like approximation (red dashed line) yields com- 
parable results for small U but overestimates the effect 
of staggering for larger values of U. 

A topological phase at half filling occurs only if we 
break the axial symmetry in the system by considering 
7 > 0. We present the non-interacting Aa^-7 phase dia- 
gram in Fig. 3(c) [7]. The interacting Xx-U phase di- 
agram for different values of 7 is shown in Fig. 3(d). 
Both semimetal and QSH phases are robust up to in- 
teractions of order U c:^ 3 — 5, at which point larger in- 
teractions drive the system into a magnetically ordered 
state. Prominently, we observe an interaction-driven NI 
to QSH transition for 7 = 0.25 and A^^ > 1.5, and a 
metal-QSH transition for 0.22 < 7 < 0.25 and A^, > 1. 

Using RDMFT for a cylinder geometry, we are able to 
directly observe the behavior of the edge states in the 
interacting system. Gapless edge states are key to differ- 
ent detection schemes of topological phases in cold-atom 
systems [23, 24]. Since topological phases are uniquely 
characterized by their helical edge states [25], a probe 
of these states is the most direct measurement [7, 26]. 
In Fig. 4, we give an example of the spectral function 
A{ky^uj) for the interaction driven NI-QSH transition at 




FIG. 4. (color online). Spectral function A(ky,uj) of interact- 
ing system clearly distinguishing between (a) NI phase with 
no edge states traversing the bulk gap at [/ = 0.5 and (b) QSH 
phase at [/ = 2 with a single pair of edge modes (per edge) 
connecting the two bulk bands. Both plots are for a = 1/6, 
7 0.25 and A^, 1.5. 

7 = 0.25, Xx = 1.5. For U = 0.5, we find no gapless 
edge modes that are connecting the two bulk bands, cor- 
responding to NI, while at /7 = 2 we clearly find a single 
pair of helical edge modes traversing the bulk gap, which 
corresponds to the QSH phase. 

Conclusion. — We have investigated the TRI 
Hofstadter-Hubbard model using RDMFT comple- 
mented by analytical arguments. We quantitatively 
determine the interacting phase diagram including two 
additional terms available in the cold- atom experiment, 
a lattice staggering and Rashba-type spin-orbit coupling. 
Interactions drive various phase transitions. Similar to 
graphene, we find that a semi- metal at half- filling turns 
into a magnetic insulator at a critical finite interaction 
strength. Rashba-type spin-orbit interactions lead to 
tunable magnetic order with collinear and spiral phases. 
We explicitly demonstrate the stability of the topolog- 
ical phases with respect to interactions, and verify the 
existence of robust helical edge states in the strongly 
correlated TI phase, which is crucial for experimental 
detection schemes. 

Note added - After the submission of this work a num- 
ber of other studies appeared [27] that investigate the 
strong-coupling limit of fermions (and bosons) in the 
presence of non-Abelian gauge fields using a spin Hamil- 
tonian similar to our Eq. (4). 
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EXPERIMENTAL DETAILS 
Realizations 

We investigate a model Hamiltonian implementing a 
non-Abelian artificial gauge field in a fermionic system. 
This model system is of general interest to topological 
insulator experiments, and directly applicable to cold-gas 
systems in optical lattices, of which a wealth of proposals 
have been made to implement these artificial gauge fields 

[!]• 

Cold-gas experiments can exhibit effective Hamiltoni- 
ans with artificial gauge fields through two main meth- 
ods: geometric phases or rotation of the cold-gas system. 
A geometry phase is created when a spatially depen- 
dent coupling connects states of an atom such that the 
Aharonov-Bohn effect is observed in the adiabatic pas- 
sage of the atom in one particular dressed state through 
a closed loop [1]. The first such realization of a geo- 
metrically induced artificial gauge field was performed 
within a continuous cold-gas of bosonic ^^Rb with a pair 
of Raman lasers and a spatially varying magnetic field 
[4]. Alternatively, the Hamiltonian of a rotation system 
can be expressed in the rotating frame, and one finds 
it is time-independent when the system is isotropic and 
introduces Coriolis terms that take the form of an artifi- 
cial gauge [2] . Such rotating systems have been achieved 
in the laboratory, with features such as a vortex lattice 
structure [3] appearing in bosonic systems. 

In a lattice, one can also introduce additional couplings 
in order to create a continuous flux and the so-called "op- 
tical flux lattice" has been proposed [5] to realize topo- 
logical non-trivial systems. However, it is not necessary 
to realize an artificial gauge over the entire lattice. In- 
stead, one can make use of the Peierls' substitution [6], 
which results in a tight-binding model that differs only by 
a phase factor in either the hopping matrix elements or, 
equivalently, of the local Wannier states. This phase fac- 
tor, the Peierls' phase, can be realized by "freezing out" 
the hopping in a lattice and then reconnecting the lattice 
sites by an additional coupling as first proposed by [7]. 
Such a realization of an Abelian artificial gauge field for a 
square lattice geometry has already been achieved in [8] 
by coupling columns of ^^Rb in a staggered optical lattice 
with Raman couplings, where the artificial magnetic flux 
alternates in sign between neighboring plaquettes. An 
extension of such a lattice to simulate a net magnetic- 
flux, which exhibits quantum Hall phases, is proposed in 
[9] via the inclusion of a super-lattice potential. Simi- 
larly, [10, 11] have shown that it is possible to generate 
non-Abelian gauge fields by asymmetric lattice shaking 
and the extension of any Abelian gauge field in a cold-gas 
experiment to a non-Abelian gauge field is possible by the 
use of an n-state atomic system [12]. The specific pro- 
posal which we follow, Goldman et al. [13], makes use of 



this idea by manipulating four states of ^Li to implement 
a non-Abelian gauge field such that a SU(2) time-reversal 
invariant form of the Hofstadter model is realized. With 
this setup, it is possible to realize any rational flux ra- 
tio a = p/q < 1/2. For large a, the edge states asso- 
ciated with these have been reported [13] to be robust 
for deviations in the flux of up to Aa ^0.01. In addi- 
tion, the non-Abelian nature of the gauge field realizes a 
spin-orbit coupling, which is the 7 term in our Hamilto- 
nian, introduced in equation (3). This term has a similar 
effect as the Rashba spin-orbit coupling in the related 
Kane-Mele system [14], a model that is closely related to 
graphene. Many other proposals exist for spin-orbit cou- 
pling in continuous ultracold gases [15-17] which could 
also be extended to lattice geometries. 

Detection and Trapping 

The major differences between cold-gas and solid-state 
experiments are the mechanisms that are available to ob- 
serve the topological properties and edge-states of the 
system, as well as the size and structure of the edge that 
is present. In optical lattices, it is not possible to per- 
form transport measurements for the Hall coefficients, 
and so one must turn to other options to directly observe 
these states. Several proposals have been made that use 
techniques such as Bragg or Raman spectroscopy [18-20] , 
which directly detect the edge states and can be applied 
regardless of the trapping potential, time-of-fiight mea- 
surements [21] to determine bulk properties, and direct 
band mapping by the use of spin-injection spectroscopy 
[17]. 

We have investigated detection using Bragg spec- 
troscopy in [18] for the non-interacting version of the sys- 
tem considered in this paper. While it is in principle pos- 
sible for us to determine the Bragg response via RDMFT 
results, we have not observed a qualitative change of the 
single-particle spectra corresponding to the edge states. 
Hence, we do not expect a qualitative deviation in Bragg 
response for excitations of these edge states compared to 
the non-interacting system. 

HERBUT'S ARGUMENT 

In Fig. 1(a) in the main text we presented numer- 
ical RDMFT results for the critical onsite interaction 
strength Uc{l/q) which marks a zero temperature quan- 
tum phase transition between a semi-metallic phase and 
a magnetically ordered phase. The system develops a 
single-particle gap inside the magnetic phase. The mag- 
netic order is antiferromagnetic and of Neel type. 

In this section, we expand our analytical analysis of the 
critical value of Uc which we obtain from a slightly modi- 
fied version of Herbut's argument. In Ref. [22] I. Herbut 
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considers the low energy field theory of graphene in the 
presence of onsite interaction U. Specifically, graphene 
exhibits four Dirac cones at low energies, a factor of two 
stemming from the two valleys (i^, K') and a factor of 
two stemming from the electronic spin. In his renormal- 
ization group analysis, Herbut extends the electronic spin 
from two to Ng fiavors and determines the /3-function of 
the Hubbard interaction to leading order in 1/Ns. In his 
notation, it reads 



(3a = -~ga-Ca~9i^O{l/Ns), 



(5) 



where Qa = Qa^, 9a = -U/{87rvFt) with hop- 

ping amplitude t and dimensionless Fermi velocity vp = 
vp/at. Here, a is the short distance cutoff of the theory 
that is set such that the Fermi velocity vf = VFdt = 1. 
The constant in the /3-function reads Ca = 2Ny, where 
Ny is the number of different valleys. For graphene, one 
finds that Ny = 2, but in the spinful Hofstadter system 
one rather finds that Ny = q and so depends on the mag- 
netic fiux per unit cell a = p/q. The critical value of the 
interaction Uc is obtained from the condition that the 
/3-function changes sign, which yields 



(6) 



2qN, 



In terms of the microscopic parameters this reads 



(7) 



where we have inserted the physical value of Ns = 2, 
in our Hofstadter model. The Fermi velocity VF{q) is 
exactly known from the non-interacting band structure 
of the Hofstadter model, and we find that it scales as 
VF{q) ^ 1/q for a = 1/q in the considered range 2 < q < 
10 [see inset in Fig. 1(a)]. As a result, we predict that 



q- 



(8) 



which agrees very well with the numerical RDMFT re- 
sults as we show in the inset of Fig. 1(a). 



DERIVATION OF EFFECTIVE SPIN 
HAMILTONIAN 

In the main text we introduce the fermionic tight- 
binding Hamiltonian including interactions and in the 
presence of the two gauge fields along the x and y links, 
which reads H = Hq -\- Hj with 



3 



h.c.} 

(9) 
(10) 



The summation j = (x, y) with x, G N runs over all 
lattice sites of the square lattice with lattice constant set 
to one. The field operator = (c^ ^, ^) is a spinor. 

At large interaction, the zeroth order Hamiltonian is 
given by the interaction part Hj. To pursue perturbation 
theory in the hopping amplitudes tx^y and derive an effec- 
tive spin Hamiltonian at half-filling up to 0{t1y/U), we 
follow the analysis in Ref. 23 and partition the complete 
Fock space into states with singly occupied sites 

S = {|ni^,ni;,n2t, . . .) : : n^^ + n^; < l} (11) 
and states with at least one doubly occupied site 

D { |ni^, ni;, n2t, . . .) : 3i : n^^ + n^; = 2} . (12) 
We partition the Hamiltonian as 



H 



A B 
C D 



PsHPs PsHPd\ 
PdHPs PdHPd) ' 



(13) 



where Ps{Pd) projects onto the subspace S{D). The 
effective Hamiltonian in the low-energy subspace S is 
obtained by projecting the resolvent operator G{E) = 
{E - H)-^ on this subspace PsG{E)Ps = [E - n{E)]-^ 
with 1-L{E) = A + B ^^jj C. Expanding to lowest or- 
der in tx^y/U and E /U yields the (energy independent) 
low-energy effective Hamiltonian 

n = PsHoPs + PsHoPd {-JjYl njrrij^)PDHoPs • 

(14) 

At half-filling the first term PsHqPs vanishes, because 
starting from subspace S each hopping event of Hq cre- 
ates a doubly occupied site and the projection onto S 
yields zero. It requires at least two hopping events for 
the system to return to the singly-occupied subspace 
5*. The second term in Eq. (14) contains all possible 
second-order virtual hopping events. Inserting the hop- 
ping Hamiltonian Hq in Eq. (14) and defining the 2x2 
complex hopping matrices = tx ex.-p{—i27r'ya^) and 
Ty{x) = ty ex.p{i27raxa^), the effective Hamiltonian can 
be written as 



n 



j,k u,fie{±x,±y} 



Ps- 



(15) 



Here, x = (1,0), y = (0, 1) are unit vectors and 

(^, M)j,fe = c]^^T^Cjcl^^T^Ck . (16) 

Note that the hopping matrices fulfill T_j, = T*. 
To arrive at Eq. (15) we have also employed that 
Pd Xlj ^j,t^j,iPD = 1, since exactly one site is doubly 
occupied in the system after the first hopping event. 

To proceed, we use that the system must return to the 
singly occupied sector 5* after the second hopping event. 
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Otherwise, the state is annihilated by the projection op- 
erator P5. It fohows that only terms with v = —fi in the 
sum in Eq. (15) are non-zero. For v = —fi = ±x, this 
demands that k = j m Eq. (16). Explicitly, we find 
(up to constants) 

(x, -x) =tlY^ cos(27r7) - icrl^, sin(27r7)] 

X Cj^rj'c\^^[5TT' cos(27r7) + icr^^, sin(27r7)]cj+i,T-/ 
= -2tl ^{cos(47r7)(5jV*5j' + 



(17) 



where the spin indices cf^g'^t^t' G {t^i} ^ind we have 
defined the spin operators Sj = ^ X^rr' ^],T^rT'Cj,T' with 
cr = {a^ ^a^ ,a^). We have used the relation 

^],c7^k,acl^^Cj^r = -"^Sj ■ Sk - ^TljUk + Uj , (18) 



where rij = 1 in our case. Since Eq. (17) is purely real, its 
hermitian conjugate term {—x^x) is identical to it. Sim- 
ilarly, we obtain for virtual hopping along the ^/-direction 

(2, -2) = -2tl 5]{cos(47rax)(5jVy^i + S^+yS]) 

3 

(19) 

where j = ix^y). The complete low-energy spin Hamil- 
tonian thus reads 



■sin(47r7) [^1^]^* " ^J^l+x] } 



ly ^< cos(47rax) 



- sin(47rax) 



(20) 



of weak, intermediate and strong interactions. RDMFT 
operates on a finite realization of a lattice by the Dyson 
equation: 



(21) 



where i and j refer to single sites of the system, is the 
non-interacting Green's function, and we adopt a vector 
notation, such that each site is decomposed into its spin- 
dependent components: 



(22) 



and similarly with the other Green's functions and self- 
energy. The core assumption of RDMFT, as with stan- 
dard DMFT, is that the self-energy is local = S^^^j. 

In the RDMFT procedure, we map each site onto 
an impurity problem by integrating out all other de- 
grees of freedom in the lattice. To solve the impurity 
problem we use a combination of exact diagonalization 
(ED), numerical renormalization group (NRG) [25] and 
the continuous-time auxiliary spin solver (CT-AUX) [26]. 
The NRG solver operates directly at T = and in real- 
frequency, which we use for the cases where 7 = 0. The 
ED and CT-AUX solvers operate in Matsubara frequen- 
cies and at finite temperature. These solvers are well 
known in solving the system without the spin-mixing due 
to the 7 term in Hamiltonian (1), and we make some 
comments about the inclusion of spin-mixing below. 

To directly observe edge states in our system, we must 
take a finite system and enforce cylindrical boundary con- 
ditions. That is, we set periodic boundary conditions in 
the y direction and open (i.e. fixed) boundary conditions 
in the x direction. In this case, we can assume ky is a 
good quantum number when no symmetry breaking has 
occurred. We investigate system sizes up to 48 sites in 
the x-direction and have allowed for symmetry breaking 
with a periodicity of up to 24 sites in the ^/-direction. 
Because we have observed no other symmetry breaking 
than AF order in the ^/-direction, we are able to restrict 
the number of impurity solver calculations to that of the 
lattice size in the x-direction. We also investigate the pe- 
riodic system to better determine the onset of magnetism 
without boundary effects. 



which is the result given in Eq. (4) of the main text. 



NUMERICAL METHOD DETAILS 



Solver details 



RDMFT details 



We have used the real-space dynamical mean-field the- 
ory (RDMFT) [24] to investigate numerically the effects 



To extend the standard Anderson impurity model to 
include spin-orbit coupling we must include, at minimum, 
a term that couples the bath orbitals of opposite spins to 
the impurity. The complete Hamiltonian for the impurity 
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is then given by 

HaIM = - ^ /i4c^ + ^ alalia 



EM. 

la 



(23) 



where Ccr is an impurity annihilation operator, a^^r is a 
bath annihilation operator, and a represents the opposite 
spin to cr. Here, the term involving Wicr represents the 
extension of the impurity model. Due to this term, the 
Weiss Green's functions for the impurity model become: 



\Vi? 



\Wi? 



Qa 



(24) 
(25) 



Formally, the CT-AUX solver is unchanged by the in- 
clusion of the spin-mixing, however one must now also 
track the spin- mixing Weiss Green's functions, Qaai for 
each configuration of auxiliary spins. This means that 
the fast-matrix updates for the solver [26] become rank- 
2 updates. In addition, we were unable to prove that 
the determinants which must be evaluated will always be 
real, and so we must consider the possibility of a "sign 
problem" (or more precisely, a phase problem) in the 
Monte-Carlo sampling. However, we have carefully an- 
alyzed all of our results and observed purely real and 
positive weights for all configurations that we consider. 



Details about gaps and magnetization 

From our converged RDMFT solutions, we determine 
several quantities directly from the output of the impu- 
rity solvers, including magnetization, double occupancy 
and the density of states at the Fermi edge for each site 
of the lattice. To identify the phase, we observe the pres- 
ence of a charge gap via several methods: a zero density 
of states at the Fermi edge, determined from extrapola- 
tion of the Matsubara Green's functions to = 0, an 
exponential decay of the imaginary time Green's func- 
tions, which is fitted to the form G(r) ~ exp(— rAgp), 
and the analytical continuation [27] to a real-frequency 
spectrum, by first continuing the self-energy [28]. 

In each case, we extract trends for different system size 
and temperature to determine the properties of the in- 
finite system at low temperatures. However, we do not 
wish to claim exact results for T = data, as this is 
beyond the capabilities of our solver, and therefore ex- 
plicitly show the results for /3 = 20. We have compared 



these results to that of ^5 = 50 for some cases and do not 
notice significant changes in the critical interactions for 
magnetic transitions. 

In our system we observe two different phase transi- 
tions in the presence of interactions: a) Semi-metal to 
magnetic order. In this transition, we observe a simulta- 
neous opening of the single-particle gap and magnetic or- 
der throughout the entire lattice, and b) Quantum spin- 
Hall phase to magnetic order. Here the system always 
possesses a single-particle gap in the bulk, and so the 
only noticeable change in the bulk is the onset of mag- 
netic order. Near the edge, the topological edge states 
also become gapped in this transition, however at larger 
temperatures we observe the single-particle gap opening 
later than the appearance of magnetization. As the mag- 
netization has destroyed time-reversal symmetry, these 
states will be susceptible to disorder and localize. The 
fate of these states in ultracold atom experiments is un- 
known and requires further investigation. 
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